/* ======Trend Inflation and the Nature of Structural Breaks in the New Keynesian Phillips Curve ======

This program estimates the unobserved components model for the NKPC with stochastic trend inflation
with no structural breaks. The variability of shocks to trend inflation are constrained to be constant
across regimes.

Written by: Pym Manopimoke

Model:

Pi_t = Pi_t^* + k Sum_{j=0}^{inf} beta^j E_{t-1}(x_{t+j}) + z_t
Pi_t^* = Pi_{t-1}^* + v_t, v_t ~iid.N(0,sig^2_v)
x_t = phi_1 x_{t-1} + phi_2 x_{t-2} + u_t, u_t|D_t~iid.N(0,sig^2_{u,D_t})
z_t = psi z_{t-1} + epsilon_t, epsilon_t ~iid.N(0,sig^2)

u_t and epsilon_t are correlated with correlation coefficient rho.

sig^2_{u,D_t} = sig^2_{u,0}(1-D_t) + sig^2_{u,1}D_t
D_t = {0 if t<= 1984Q3, 1 otherwise}

========================================================================================================*/

new;
library optmum,PGRAPH;
pqgwin auto;
cls;
 
load data[250,11]=data.prn;  
data=data[2:250,.];
t=rows(data);
@Set the sample to 1952:1-2007:3@
startsmpl=21;
endsmpl =243; 

qtr= data[startsmpl:endsmpl,1];  
/*Column 2 for GDP chained-weighted inflation, 7 for CPI inflation, 8 for PCE inflation*/  
inf=data[startsmpl:endsmpl,8];  @PCE inflation@
inff1=data[startsmpl+1:endsmpl+1,8];  @one-quarter ahead inflation@

xx=data[startsmpl:endsmpl,5];  @CBO output gap@
t=rows(inf~xx);

@Define dummy variable series D_t to account for the Great Moderation@
GM = 130; @1984:2@
dmat = ones(GM,1)|zeros(t-GM,1); 


@================= Initialize Global Variables============@

tvp_no=4; 
start= 11;
t1=75;t2=125;
  
prmtr_in= 1.0297822 |
          0.45991641 |
          0.34008988 |
          0.73952899 |
          -0.083434780 |
          1.4732290 |
          -4.3488046 |
          1|
          0.018851190 ;
     
{xout,fout,gout,cout}=optmum(&lik_fcn,PRMTR_IN);

PRM_FNL=TRANS(xout);             @ Estimated  coefficients, constrained@
xout;
hout0=hessp(&lik_fcn,xout);
hout0;
hout=inv(hout0);
grdn_fnl=gradfd(&TRANS,xout);
Hsn_fnl=grdn_fnl*hout*grdn_fnl';
SD_fnl =sqrt(diag(Hsn_fnl));     @Standard errors of the estimated coefficients@
 
{B_TT,stdprederr} =FILTER(XOUT);
_pdate=0;
_pltype = {4 6 3 3};
title("Inference on inflation trend");
_plegctl={2 3.8 5.7 5};
_plegstr="Actual Inflation\000Trend Inflation\00090\% Confidence Bands";
xy(qtr[start:t],inf[start:t]~b_tt[.,1] );
 
output file=UC1.out reset;

"== FINAL OUTPUT ===============================";
"Likelihood value is ";; -fout;
"Estimated parameters (s.e.) tstat are:"; 
prm_fnl~sd_fnl~prm_fnl./sd_fnl;
"===============================================";
output off;
end;

@ end OF MAIN PROGRAM @
@========================================================================@

PROC LIK_FCN(PRMTR1);

local PRMTR,p00,p11, p22,siguGM0,siguGM1,sigu,sigv,sigep,rho,psi,PHI1,PHI2,k,
xfor,dmt,ome,ome0,ome1,ome2,subf,subq,H,vv,vecp,FF,A,yt,
PROB__0,PROB__1,PROB__2,LIKV,J_ITER,Q0,Q1,Q2,F0,F1,F2, 
B_TL00,B_TL01,B_TL02,B_TL10,B_TL11,B_TL12,B_TL21,B_TL20,B_TL22,     
P_TL00,P_TL01,P_TL02,P_TL10,P_TL11,P_TL12,P_TL20,P_TL21, P_TL22, 
F_CAST00,F_CAST01,F_CAST02,F_CAST10,F_CAST11,F_CAST12,F_CAST20,F_CAST21,F_CAST22,
SS00,SS01,SS02,SS10,SS11,SS12,SS20,SS21,SS22,
B_LL0,B_LL1,B_LL2,P_LL0,P_LL1,P_LL2,pr_ll0_,pr_ll1_pr_ll2_, 
B_TT00,B_TT01,B_TT02,B_TT10,B_TT11,B_TT12,B_TT20,B_TT21,B_TT22,
P_TT00,P_TT01,P_TT02,P_TT10,P_TT11,P_TT12,P_TT20,P_TT21,P_TT22,
PR_VL00,PR_VL01,PR_VL02,PR_VL10,PR_VL11,PR_VL12,PR_VL20,PR_VL21,PR_VL22,
PR_VAL,PRO_00,PRO_01,PRO_02,PRO_10,PRO_11,PRO_12,PRO_20,PRO_21,PRO_22;  


EXTERNAL PROC TRANS, V_PROB;
 
           PRMTR=TRANS(PRMTR1);
           siguGM0 = prmtr[1,1];
           siguGM1 = prmtr[2,1];
           sigv= prmtr[3,1];  
           sigep = prmtr[4,1];  
           rho = prmtr[5,1];
           phi1 = prmtr[6,1];
           phi2 = prmtr[7,1];
           psi=prmtr[8,1];    
           k=prmtr[9,1];

           A=1~0;
           H = 1~1~0~0|0~0~1~0;
           
           B_LL0 = inf[1]|0|xx[1]|0;      @initial guess for coefficients@
           P_LL0 = eye(4)*10000;      @uncertainty associated with initial guess@
           
           subF = phi1~phi2| 1~0; 
           subQ= siguGM0^2~0|0~0;
           vecp = inv(eye(4)-subF.*.subf)*vec(subq);
           vv = reshape(vecp',2,2);
           p_ll0[3:4,3:4] = vv;
           p_ll0[2,2] = sigep[1,1]^2/(1-psi^2);
         
/*  START ITERATION */
/*===================================================================== */

       LIKV=0.0;
       J_ITER = 1;
       do until J_ITER>T;

            yt=inf[j_iter]|xx[j_iter];
     
            dmt = dmat[j_iter];
            sigu=siguGM0*dmt + siguGM1*(1-dmt);
     
            Q0= sigv^2~0~0~0|
                0~sigep^2~rho*sigu*sigep[1,1]~0|
                0~rho*sigu*sigep~sigu^2~0|
                0~0~0~0;
           F0 = 1~0~0~0 |
                0~psi~0~0|
                0~0~phi1~phi2|
                0~0 ~1~0;
           FF= phi1~phi2|1~0;
           xfor = k*A*(inv(eye(2)-FF));   
           ome0 =0~0~xfor| 0~0~0~0;
           
           
        /* KALMAN FILTER */
      
           B_TL00 = F0 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=0@
           P_TL00 = F0* P_LL0 * F0' + Q0;
           F_CAST00= yt - H * B_TL00 - ome0* B_TL00;
           SS00= H * P_TL00 * H' ; @VARIANCE OF FORECAST ERROR@
           B_TT00 = B_TL00 + (P_TL00 * H' *inv(SS00)) * F_CAST00;
           P_TT00 = (EYE(tvp_no) - (P_TL00 * H' *inv(SS00) * H )) * P_TL00;
            
         
       /* HAMILTON FILTER */ 
    
           PR_VL00=V_PROB(F_CAST00,SS00) ;    @Pr[St,Yt/Yt-1]@
           PR_VAL=PR_VL00 ;
           
       /* COLLAPSING TERMS */
      
           B_LL0= B_TT00;
           P_LL0= (P_TT00+(B_LL0-B_TT00)*(B_LL0-B_TT00)' );
           if J_ITER < START; goto skip; endif;
           LIKV = LIKV -LN(PR_VAL);
           skip:
       J_ITER = J_ITER+1;
       ENDO;

/*  END ITERATION    */


RETP(LIKV);
ENDP;

@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
proc TRANS(prmtr1); @ constraining values of reg. coeff.@
local prm_new,d0,d1;
prm_new=prmtr1; 
 
     prm_new[1:4,1] = sqrt(prmtr1[1:4,1]^2);
     prm_new[5] =prmtr1[5,1]./(1+abs(prmtr1[5,1]));
     d0 =prmtr1[6,1]./(1+abs(prmtr1[6,1]));
     d1 =(1-abs(d0))*prmtr1[7,1]./(1+abs(prmtr1[7,1]))+ abs(d0)-d0^2;
     prm_new[6,1]=2*d0;
     prm_new[7,1]=-1*(d0^2 + d1) ;
     prm_new[8,1] =prmtr1[8,1]./(1+abs(prmtr1[8,1]));
    
retp(prm_new);
endp;


@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
PROC V_PROB(EV, HE);      @ CALCULATES    Pr[Yt/St,Yt-1]@

local VAL;
VAL=(1/SQRT(2*PI*det(HE)))*EXP(-0.5*EV'inv(HE)*EV);

RETP(VAL);
ENDP;
@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@


proc (2) = FILTER(PRMTR1);

local PRMTR,p00,p11, p22,siguGM0,siguGM1,sigu,sigv,sigep,rho,psi,PHI1,PHI2,k,
xfor0,xfor1,xfor2,dmt,ome,ome0,ome1,ome2,subf,subq,H,vv,vecp,FF,A,yt,confint,lowerCI,upperCI,stderr,stdprederr,xforecast,
PROB__0,PROB__1,PROB__2,LIKV,J_ITER,Q0,Q1,Q2,F0,F1,F2,B_M,PR_TTM,PR_TLM,twostep,
pr_tl00,pr_tl01, pr_tl02, pr_tl10,pr_tl11,pr_tl12,pr_tl20,pr_tl21,pr_tl22,
B_TL00,B_TL01,B_TL02,B_TL10,B_TL11,B_TL12,B_TL21,B_TL20,B_TL22,     
P_TL00,P_TL01,P_TL02,P_TL10,P_TL11,P_TL12,P_TL20,P_TL21, P_TL22, 
F_CAST00,F_CAST01,F_CAST02,F_CAST10,F_CAST11,F_CAST12,F_CAST20,F_CAST21,F_CAST22,
SS00,SS01,SS02,SS10,SS11,SS12,SS20,SS21,SS22,
B_LL0,B_LL1,B_LL2,P_LL0,P_LL1,P_LL2,pr_ll0_,pr_ll1_pr_ll2_, 
B_TT00,B_TT01,B_TT02,B_TT10,B_TT11,B_TT12,B_TT20,B_TT21,B_TT22,
P_TT00,P_TT01,P_TT02,P_TT10,P_TT11,P_TT12,P_TT20,P_TT21,P_TT22,
PR_VL00,PR_VL01,PR_VL02,PR_VL10,PR_VL11,PR_VL12,PR_VL20,PR_VL21,PR_VL22,
PR_VAL,PRO_00,PRO_01,PRO_02,PRO_10,PRO_11,PRO_12,PRO_20,PRO_21,PRO_22,
pi2hat_00,pi2hat_01,pi2hat_02,pi2hat_10,pi2hat_11,pi2hat_12,
pi2hat_20,pi2hat_21,pi2hat_22,B_T2L00,B_T2L01,B_T2L02,
B_T2L10,B_T2L11,B_T2L12,B_T2L20,B_T2L21,B_T2L22;  

EXTERNAL PROC TRANS, V_PROB;


           PRMTR=TRANS(PRMTR1);
 
           siguGM0 = prmtr[1,1];
           siguGM1 = prmtr[2,1];
           sigv= prmtr[3,1];  
           sigep = prmtr[4,1];  
           rho = prmtr[5,1];
           phi1 = prmtr[6,1];
           phi2 = prmtr[7,1];
           psi=prmtr[8,1];    
           k=prmtr[9,1];

           A=1~0;
           H = 1~1~0~0|0~0~1~0;
           
           B_LL0 = inf[1]|0|xx[1]|0;      @initial guess for coefficients@
           P_LL0 = eye(4)*10000;      @uncertainty associated with initial guess@
           
           subF = phi1~phi2| 1~0; 
           subQ= siguGM0^2~0|0~0;
           vecp = inv(eye(4)-subF.*.subf)*vec(subq);
           vv = reshape(vecp',2,2);
           p_ll0[3:4,3:4] = vv;
           p_ll0[2,2] = sigep[1,1]^2/(1-psi^2);

           B_M=ZEROS(T,4);  @will store (1,1)element of B_TT@
           stdprederr=zeros(t,1);
 
/*  START ITERATION */
/*===================================================================== */

       LIKV=0.0;
       J_ITER = 1;
       do until J_ITER>T;
            yt=inf[j_iter]|xx[j_iter];
     
            dmt = dmat[j_iter];
            sigu=siguGM0*dmt + siguGM1*(1-dmt);
     
            Q0= sigv^2~0~0~0|
                0~sigep^2~rho*sigu*sigep[1,1]~0|
                0~rho*sigu*sigep~sigu^2~0|
                0~0~0~0;
           F0 = 1~0~0~0 |
                0~psi~0~0|
                0~0~phi1~phi2|
                0~0 ~1~0;
           FF= phi1~phi2|1~0;
           xfor0 = k*A*(inv(eye(2)-FF));   
           ome0 =0~0~xfor0| 0~0~0~0;
             
        /* KALMAN FILTER */
      
           B_TL00 = F0 * B_LL0;  @WHEN S_{t-1}=0, S_{t}=0@
           P_TL00 = F0* P_LL0 * F0' + Q0;
           F_CAST00= yt - H * B_TL00 - ome0* B_TL00;
           SS00= H * P_TL00 * H' ; @VARIANCE OF FORECAST ERROR@
           B_TT00 = B_TL00 + (P_TL00 * H' *inv(SS00)) * F_CAST00;
           P_TT00 = (EYE(tvp_no) - (P_TL00 * H' *inv(SS00) * H )) * P_TL00;
           
           B_TT=B_TT00 ;
           B_M[J_ITER,.]=B_TT'; 
           stderr= (f_cast00[1,1]/sqrt(ss00[1,1]));
           STDPREDERR[j_iter,.] = stderr;
       /* HAMILTON FILTER */ 
    
           PR_VL00=V_PROB(F_CAST00,SS00) ;    @Pr[St,Yt/Yt-1]@
           PR_VAL=PR_VL00 ;
           
       /* COLLAPSING TERMS */
           
           B_LL0= B_TT00;
           P_LL0= (P_TT00+(B_LL0-B_TT00)*(B_LL0-B_TT00)' );
           if J_ITER < START; goto skip; endif;
           LIKV = LIKV -LN(PR_VAL);
           skip:
 
       J_ITER = J_ITER+1;
       ENDO;

/*  END ITERATION    */
/*===================*/

retp(B_M[START:T,.],stdprederr[start:t,.]);
ENDP;


@==========================================================@

